function [auxOut,f,errorMes] = numF(params)
 
% Computing the steady 
[auxOut,errorMes] = DSGEmodel_ss(params);
params = auxOut.params;
xss = auxOut.xssTrans;
yss = auxOut.yssTrans;
 
%Unfold params
BETTA= params.BETTA;
B= params.B;
CHI= params.CHI;
CHI0= params.CHI0;
THETA= params.THETA;
DELTA= params.DELTA;
ALFA= params.ALFA;
PHI= params.PHI;
PHIzero= params.PHIzero;
ZI= params.ZI;
KAPAw= params.KAPAw;
ETA= params.ETA;
RHOR= params.RHOR;
PHIpai= params.PHIpai;
PHIpai_1= params.PHIpai_1;
PHIy= params.PHIy;
PHIy_1= params.PHIy_1;
PHIc= params.PHIc;
PHIc_1= params.PHIc_1;
PHIl= params.PHIl;
NU= params.NU;
U0= params.U0;
U0d= params.U0d;
RHOz= params.RHOz;
RHOd= params.RHOd;
RHOn= params.RHOn;
RHOp= params.RHOp;
RHOa= params.RHOa;
Kss= params.Kss;
OUTPUTss= params.OUTPUTss;
PAIss= params.PAIss;
MUZss= params.MUZss;
Css= params.Css;
AA= params.AA;
lss= params.lss;
Rss= params.Rss;
Wss= params.Wss;
 
% Current values in the model;
% x;
c_ba1 = xss(1);
muz_cu = xss(2);
d_cu = xss(3);
n_cu = xss(4);
paistar_cu = xss(5);
a_cu = xss(6);
% xp;
c_ba1p = xss(1);
muz_cup = xss(2);
d_cup = xss(3);
n_cup = xss(4);
paistar_cup = xss(5);
a_cup = xss(6);
 
% y
c_cu = yss(1);
r_cu = yss(2);
pai_cu = yss(3);
w_cu = yss(4);
l_cu = yss(5);
output_cu = yss(6);
mc_cu = yss(7);
evf_cu = yss(8);
vf_cu = yss(9);
p1_cu = yss(10);
p1Real_cu = yss(11);
dc_cu = yss(12);
% yp;
c_cup = yss(1);
r_cup = yss(2);
pai_cup = yss(3);
w_cup = yss(4);
l_cup = yss(5);
output_cup = yss(6);
mc_cup = yss(7);
evf_cup = yss(8);
vf_cup = yss(9);
p1_cup = yss(10);
p1Real_cup = yss(11);
dc_cup = yss(12);
 
%% Setting the dimension of f matrices 
f  = zeros(length(xss)+length(yss),1);
 
%% f Function evaluation
% START DISPLAYING f
f(1,1)=exp(-vf_cu)*(exp(d_cu)*(((exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))/Css^CHI0)^(1 - CHI)/(CHI - 1) - U0d + (PHIzero*exp(n_cu)*(1 - exp(l_cu))^(1 - 1/PHI))/(1/PHI - 1)) - U0 + (AA*BETTA)/exp(evf_cu)^(1/(ALFA - 1))) - 1;
f(2,1)=exp(-w_cu)*(KAPAw*Wss - (Css^CHI0*PHIzero*exp(n_cu)*((exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))/Css^CHI0)^CHI*(KAPAw - 1))/(1 - exp(l_cu))^(1/PHI)) - 1;
f(3,1)=(BETTA*exp(-d_cu)*exp(-pai_cup)*exp(d_cup)*exp(r_cu)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(4,1)=- (exp(-a_cu)*exp(-mc_cu)*exp(w_cu)*exp(l_cu)^THETA)/(Kss^THETA*(THETA - 1)) - 1;
f(5,1)=1 - (exp(-mc_cu)*(ETA + (ZI*exp(pai_cu)*(exp(pai_cu)/PAIss^NU - 1))/PAIss^NU - (BETTA*ZI*exp(-d_cu)*exp(-output_cu)*exp(d_cup)*exp(muz_cup)*exp(output_cup)*exp(pai_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*(exp(pai_cup)/PAIss^NU - 1)*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(PAIss^NU*(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI) - 1))/ETA;
f(6,1)=Rss*exp(PHIl*log(exp(l_cu)/lss) + PHIc*(log(exp(muz_cu)) - log(MUZss) - exp(c_ba1) + exp(c_cu)) + PHIc_1*(log(exp(muz_cup)) - log(MUZss) - exp(c_ba1p) + exp(c_cup)) - PHIpai*(log(exp(paistar_cu)) - log(exp(pai_cu)) + log(PAIss)) - PHIpai_1*(log(exp(paistar_cup)) - log(exp(pai_cup)) + log(PAIss)) + PHIy*log(exp(output_cu)/OUTPUTss) + PHIy_1*log(exp(output_cup)/OUTPUTss)) - exp(r_cu);
f(7,1)=(exp(output_cu)*((ZI*(exp(pai_cu)/PAIss^NU - 1)^2)/2 - 1))/(exp(c_cu) + DELTA*Kss) + 1;
f(8,1)=Kss^THETA*exp(-output_cu)*exp(a_cu)*exp(l_cu)^(1 - THETA) - 1;
f(9,1)=exp(-evf_cu)*((exp(vf_cup)*exp(muz_cup)^((CHI - 1)*(CHI0 - 1)))/AA)^(1 - ALFA) - 1;
f(10,1)=(BETTA*exp(-d_cu)*exp(-p1_cu)*exp(-pai_cup)*exp(d_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(11,1)=(BETTA*exp(-d_cu)*exp(-p1Real_cu)*exp(d_cup)*exp(muz_cup)^(CHI*(CHI0 - 1) - CHI0)*(exp(c_cu) - B*exp(-muz_cu)*exp(c_ba1))^CHI*((AA*exp(-vf_cup))/(exp(evf_cu)^(1/(ALFA - 1))*exp(muz_cup)^((CHI - 1)*(CHI0 - 1))))^ALFA)/(exp(c_cup) - B*exp(-muz_cup)*exp(c_cu))^CHI - 1;
f(12,1)=log(exp(c_cu)) - log(exp(c_ba1)) - dc_cu + log(exp(muz_cu));
f(13,1)=exp(-c_ba1p)*exp(c_cu) - 1;
f(14,1)=RHOz*log(exp(muz_cu)/MUZss) - log(exp(muz_cup)/MUZss);
f(15,1)=RHOd*log(exp(d_cu)) - log(exp(d_cup));
f(16,1)=RHOn*log(exp(n_cu)) - log(exp(n_cup));
f(17,1)=RHOp*log(exp(paistar_cu)) - log(exp(paistar_cup));
f(18,1)=RHOa*log(exp(a_cu)) - log(exp(a_cup));
% END DISPLAYING f
 
end
